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This paper presents an optimal control method for a class of distributed-parameter systems governed by 
first order, quasilinear hyperbolic partial differential equations that arise in many physical systems. Such 
systems are characterized by time delays since information is transported from one state to another by wave 
propagation. A general closed-loop hyperbolic transport model is controlled by a boundary control embed- 
ded in a periodic boundary condition. The boundary control is subject to a nonlinear differential equation 
constraint that models actuator dynamics of the system. The hyperbolic equation is thus coupled with the 
ordinary differential equation via the boundary condition. Optimality of this coupled system is investigated 
using variational principles to seeft an adjoint formulation of the optimal control problem. The results are then 
applied to implement a model predictive control design for a wind tunnel to eliminate a transport delay effect 
that causes a poor Mach number regulation. 


I. Introduction 


The terminology of distributed-parameter system is used to describe physical systems whose state variables are 
functions of both space and time that are usually modeled by partial differential equations (PDEs). Optimal control 
of distributed-parameter systems has been studied extensively in mathematical literature, but applications of this type 
of control are limited in practical control engineering. Hyperbolic PDEs are used to model transport systems whose 
information is carried from one point to another within those systems as a function of space and time. 1 Examples 
of transport systems are numerous in many physical applications such as fluid flow in gas distribution pipelines, 2 air 
traffic systems, 3 highway traffic systems, 4 to name a few. These equations describe wave propagation that exists in 
transport systems to propagate information from one point to another within the continuum. A characteristic of a 
hyperbolic system is time delay which describes a finite time that information is transported from one location to 
another. Consider the simplest hyperbolic equation 


dy 

dt 



= 0 


CD 


The solution of this equation is 

y = ( 2 ) 

Equation (2) describes the solution of a wave of a transport process that propagates at a wave speed a and the term 
x/a is a variable transport time delay. 

As with any PDEs, boundary conditions are used to specify configurations of these transport systems. If the 
information is carried in one direction without returning to its starting position, then we say that the system is open- 
loop. An example of an open-loop system is gas flow through an aircraft engine. On the other hand, if the information 
returns to its starting position, then the system is said to be closed-loop. An example of a closed-loop system is the 
cardio- vascular circulatory system. Boundary conditions associated with closed-loop systems are usually periodic in 
nature. 

The flow of information is usually supplied at the system boundary by a forced process that provides a motive 
force to move the information along the way by wave propagation. For example, a common device for accomplishing 
this objective in fluid transport systems is a pump which supplies a positive pressure head to displace the fluid volume 
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in the flow direction. Such a process whereby the control is applied at the boundary of the continuum is called a 
boundary control process. In real systems, boundary control processes are often governed by other auxiliary dynamical 
processes. For example, a positive-displacement pump may be driven by an electrical motor that imposes a constraint 
on the pump speed according to the motor torque dynamics. The pump speed in this case is a boundary control 
variable. 

In this study, we present an optimal control method for a class of closed-loop hyperbolic distributed-parameter 
systems whose boundary control process is constrained by a dynamical system governed by a nonlinear ordinary 
differential equation (ODE) via a periodic boundary condition. The optimal control problem is formulated by an 
adjoint method to seek necessary conditions for optimality via calculus of variations. The theory is then used to 
develop a model predictive linear-quadratic regulator optimal control that results in a modified Riccati equation. The 
model predictive control computes control inputs based on model information at the current time that need to be apply 
ahead of error signals in order to minimize system disturbances due to time delays. Thus, it is effectively a feedforward 
control. We apply this theory to solve a flow control problem in a wind tunnel. 


II. Problem Formulation 


Transport phenomena are governed by the conservation laws equations which dictate the conservation of some 
quantities such as traffic density, mass flux, and enthalpy. These equations are generally hyperbolic in nature. For a 
1-D system, these equations are expressed in a conservation form as 5 

% + + Q (y> *) = 0 ■ Vx € [°> £ j . t € [0, tfl ' (3) 

where y (x, t) : [0, L] x [0, tf] ® n in class C 1 is a vector of conserved quantities, F (y, x) is a flux function, and 
Q (y, x) is a non-homogeneous source term. 

By explicit differentiation, Eq. (3) can be rewritten in as 

% + A (y..*) + B (y> x ) = 0 Vz e [o, L] , t e [0, t f ] (4) 

where A (y, x) : E 71 x [0, L] — > R n x R n is a characteristic matrix such that A (y, x) — F y (y , x) and B (y , x) : 
R n x [0, L] R n is a non-homogeneous source term such that B (y , x) ~ Q (y , x) + F x (y, x). 

Equation (4) is a system of first order, quasilinear hyperbolic equations due to the fact that the matrix A has n real, 
distinct eigenvalues such that 


.(A) < A 2 (A) ••<...< A m (A) < A m+ i (A) < A m 4.2 <!...< A n (A) 

for all y (x y t) € M n , x € [0, L], and t £ [0, tj], where m < n are number of negative eigenvalues of A. 

The eigenvalues are the wave speeds of the transport system and the direction of the wave propagation is called a 
characteristic direction. If m > 0, the information in the continuum is carried in both the upstream and downstream 
directions by negative wave speeds A ? , i — 1,2,. . . , m and positive wave speeds A u i ~ m -f 1, m 4- 2, . , . , n; 
respectively. If the solution domain is 0 <x < L, then for the information to be transported in the upstream direction 
by the negative wave speeds, information must exist at the boundary x — L. Similarly, information must also exist at 
the boundary x — 0 in order for the information to be carried downstream by the positive wave speeds. Therefore, the 



This is known as the boundary condition compatibility. 


In a closed-loop transport system, information is carried from one point to another and then returned back to 
the starting position as illustrated in Fig. 1. To enable this information flow, a periodic boundary control process is 
embedded within the system. For a closed-loop system, the boundary conditions at x = 0 are affected by the boundary 
conditions at x = L since the information must be returned to its starting position. Thus, in general for a closed-loop 
system, we specify the following general nonlinear periodic boundary condition for Eq. (4) 

y (0, t) = g (y (X, t ) , u {t)) 'it 6 [0, tf] (5) 

where u ( t ) : [0, tf] — > R m in class G 1 is a boundary control vector, and g (y (L, t),u) : 1" x R m — > R n is a 
nonlinear forcing function that relates the transport state vectors at x — 0 and x — L and the boundary control vector 
u. 
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Fig. 1 - Closed-Loop Transport Model 

To ensure the boundary condition compatibility for all signs of eigenvalues, the Frechet derivative or the Jacobian 
of g with respect to y (L,t) is required to have a full rank or 

dim [g y(i , t) ] = n (6) 

In many real systems, boundary control processes are actually controlled by other auxiliary processes. These 
auxiliary processes may be dynamical so that their dynamics can be described by a nonlinear ODE 

u - f (y (0, t) , y (L, t) , u, v) (7) 

where v (t) : [0, £/] — + W l is an auxiliary control vector that belongs to a convex subset of admissible auxiliary control 
Vod Q K 1 , and f (y (0, t) , y (L, £) ,u,v) : x R n x R m x R l — ► M m is a nonlinear function. 

Thus the auxiliary control vector v actually influences the boundary control vector u, which in turn controls the 
behavior of the closed-loop transport system described by Eq. (4) and the periodic boundary condition (5). 


HI. Optimal Control Theory 

Optimal control and optimization theories of hyperbolic systems has been studied extensively in mathematical 
literature. Within the theoretical framework of systems governed by PDEs, control of such systems can exist as dis- 
tributed control, boundary control, interior pointwise control, or others. Hou and Yan studied the long time behavior 
of solutions for an optimal distributed control problem for the Navier-Stokes equations. 5 Nguyen et al investigated a 
flow control problem with interior pointwise control. 7 Optimal control problems of transport systems with boundary 
control have been examined for many different types of constraints imposed on either state or control variables. Ray- 
mond and Zidani investigated necessary optimality conditions in the form of a Pontryagin’s minimum principle for 
semilinear parabolic equations with pointwise state constraints and unbounded control. 8 Casas et al established second 
order sufficient conditions for local optimality of elliptic equations with pointwise constraints on the boundary control 
and equality and set-constraints on state variables. 9 Kazemi obtained adjoint equations for a degenerate hyperbolic 
equation. 10 

Adjoint method is well-known in optimal control and optimization theories as it provides an indirect method 
for solving optimization problems. For a transport system governed by hyperbolic equations, two types of adjoint 
formulation are used: discrete adjoint and continuous adjoint Any hyperbolic equation can generally be discretized 
j into a system of ODEs by means of numerical discretization techniques such as finite-difference or finite-element 

j methods. If the adjoint method is formulated with the discretized hyperbolic equations, then this is known as discrete 

} adjoint method. On the other hand, continuous adjoint method is normally applied directly to the original hyperbolic 

equations. This method has been used in aerodynamic design optimization studies involving the Euler and Navier- 
| Stokes equations. 11 Nadarajah. et al compared the discrete and continuous adjoint methods in aerodynamic design 

I optimization and suggested that the continuous adjoint method affords a certain advantage over the discrete adjoint 

| method for Navier-Stokes flow problems . 12 In the present study, we apply this method to the present hyperbolic sy stem 

with nonlinear differential equation constraints on a periodic boundary condition. To formulate the continuous adjoint 
| method for this system, we seek a solution of the hyperbolic system above that minimizes the following multi-objective 

I cost functional ^ 

| min J(y,u,v) = f f Li(y)dxdt + [ L 2 (y° > y' L ,u,v) dt (8) 

J o Jo Jo 
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where L\ is an objective function defined over the continuum of the transport system, £2 is an objective function 
defined on the system boundary, and the superscripts 0 and L are short hand notations denoting the associated vector 
quantity evaluated at x = 0 and x = L, respectively. 

The following assumptions are required: 

( A1 ) : Eq. (4) admits smooth solutions for shock-free conditions. 

(A2) : v £ £ 2 > tihe space of real value functions in R 1 for which the norm ||v|| is square-integrable. 

(A3): The Frechet derivatives of L 2 , B, g, and f with respect to y, y°, y 1 ', u, and v exist and are bounded 
so as to satisfy the Lipschitz condition. 

We note that Eq. (4) also has discontinuous solutions known as entropy solutions i4 which will not be treated here. 

The transport system above is posed as a boundary control problem of hyperbolic equations with nonlinear differ- 
ential equation constraints. 

Lemma 1: Let D be a nonlinear differential operator and D* be its adjoint differential operator such that for some 
z (x, t) G M n and A (x, t) G E n 


n a dz 
Dz = A— 

ox 

D -a = |(a t x) 


where the superscript T is the transpose operator, then the following inner product operation in £2 is equivalent 


(Dz,\} {x>t) = ~ (z,D* A) {S;t) + {z L , (A T A) L } t - (z°, (A t A)‘ 


(9) 


Proof: The inner product ( Dz , A) in £2 is 


Integrating by parts yields 


( Dz > X )(x,t) 



\ T Dzdxdt 



\ r Dzdxdt = 


[ f [ z T D*\dxdt + f ' [(z L ) T (A L ) T \ L - (z°) T (A°) T A 0 ] dt 
Jo Jo Jo ^ -* 


We define the following inner products 


(A - (z», (A^)°), = I’’ [^) T (A‘) T A* - (,") T (A«) T A«] * 

Equation (9) is thus obtained. 

Definition 1: Let F : X — ► Y be a functional with X, Y in Banach spaces and a G X. If there exists a continuous 
linear operator VF (a) : X F for any variation 5 G X such that 


lim 

e— >0 


VF (a) S ~ 


Ffat + eS) -F(g} 
e 


0 


then VF (a) is called a Gateaux derivative of F at a. 16 
Definition 2: The following Hamiltonians are defined 


Hi (y, x, A) = Lj - A t B (10) 

Hi (y°, y L , u, v,p) = L 2 + pt r £ (1 1) 

We are now ready to state the necessary conditions for optimality. 

Theorem: If (A1)-(A3) are fulfilled and if (y, u,v) is an optimal solution of Eq. (8), then there exist adjoint 
variables A (x, t) : [0, L] x [0, tf] — ► E n and /x : [0, tf] — *■ R k that satisfy the following dual adjoint system 


(A T A) i 


A t + (A T A) a 


H ly 


L + g^L 


'rHZ y = 0 

(A'X)° + gJ L Hj ty0 


A (Xytf) = 0 
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A - -Hi u - gj (A ! A)° - slHjy, 

(x (t f ) = 0 

with a terminal time transversality condition 

J q L Li\t=t f dx+ L 2 \ t=t/ =0 (14) 

such that the optimal control is one that satisfies the foUo wing Pontry agin ’s minimum principle 

v = arg mill H 2 (y°,y L , u, v, (j) (15) 

VSVad 

Proof: Let a = (z, p) be solutions to /3 = (y, u) in variations for a variation q in v, then the variation in the cost 
functional from Eq. (8) is computed as 

A J (q) = V J 08) + J 08, v + q) - J 08, v) > 0 

where the Gateaux derivative of J at j3 is evaluated as 



W(/3) * (Hl y ,z) {x>t) - (A, z t ) (jSit) - (A,^) (Xit) + (^ y o lZ °) t 

+ {^y L}2,L \ + ~ + 

From the boundary condition (5), we have the following variations 

Z 0 = g^Z 1, +g^U 

From Lemma 1 and the variations in the boundary condition (5) plus vanishing variations in initial conditions for 
Eqs. (4) and (7), this becomes 


VJ (0) = <A t + D*A + Hj >y , (x, i,) , A (x,tf)) x +(g^Hl y o + g; L (A 1 A)° + 23^ - (A t A) £ , z) t 

+ (gj (A T A)° + gJfr ? Ty° + + A t , p) t - m t (if) p (tf) +st ^ + r 2 | t=: 

Setting VJ (/3) = 0 for arbitrary variation a results in Eqs. (12)-(14). Then the variation in the cost functional 
becomes 

AJ(q )=/ [H 2 (y°,y L , u,v + q,n) -/x T u]dt- I [H 2 (y° ,y L ,u,v, fx) - fx T ti] dt > 0 

Jo Jo 

This leads to the Pontry agin J s minimum principle 

H 2 (y°,y £ ,u,v + q ,n) > H 2 (y°,y L ,u,v,p) 


for all values of the variation q. 

Equation (15) is the equivalent Weierstrass condition for strong variations. For weak variations when v is uncon- 
strained and L\ and £3 are convex, the Pontry agin’ s principle leads to the Legendre -Clebsch condition 15 


H 2l v (y°,y L ,u,v,0) = 0 1 
H 2>vv (y°, y L , u, v, fx) > 0 J 
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IV. Linear-Quadratic Optimal Control for Scalar Hyperbolic Systems 

Linear advection hyperbolic equations are used to model some transport processes that involve small- amplitude 
wave propagation in only one direction such as pressure wave propagation in a fluid continuum. Our objective is to 
develop a linear-quadratic regulator (LQR) for the following linear advection equation 

+ S + + = 0 (17) 

where y t) El is a transport state variable, a (x) > 0 is the wave speed, /3 (x> t) is a dissipative term, and w (a:, t) 
is a disturbance. 

Equation (17) is subject to a zero initial condition and a periodic boundary condition for a closed-loop system 


y(0 > t)*=.G.u(t) + ITy{L > t) (18) 

where u (t) £ E m is a boundary control vector, G:M->lx M m is a constant-valued matrix, and H is a constant. 
The following linear dynamics is imposed on the boundary control vector u with a zero initial condition 


u - Cu 4 Dv 4- E y (0, t) 4 'Ey (L, t) (19) 

where C : R — > 3R m xl m is a constant- valued state transition matrix, D : M ■— > 3R m x R l is a constant-valued control 
transition matrix, and E and F are constants. 

We want to minimize the following linear-quadratic cost functional with respect to v for a fixed final time if 

( 20 ) 

( 21 ) 

( 22 ) 

(23) 

with the transversality conditions A (x,tf) '= 0 and jx [if) ~ 0 along with the stationary condition obtained from Eq. 
(16) 


min J 


r*s 

Jo 


\py 2 (0, t) + ~u T Qu + |V T Rv 


di 


where P > 0, Q > 0, and R > 0 are weighting factors. 

The dual adjoint systems from (12) and (13) for the optimal control problem are given as 

■f- 

a (L) X ( L,t ) = Ha (0) X (0, t) + HPy (0, t) + (F t + HB t ) /j, 
it = -Qu - (C T + G t E t ) jx - G T a (0) A (0, t) - G T Py (0, t) 


Rv +G T (j, — 0 (24) 

Equations (!?)-( 19) and (21)-(23) form a two-point boundary value PDE-ODE problem. Even though the PDEs 
are linear and scalar, the two-point boundary conditions ( 1 8) and (22) pose a challenge in obtaining a general feedback 
control solution. To see this, we solve Eq. (17) using the characteristic method which yields the following solution 


y (x, t ) = 


o t < u {x) 

a (x, t) [f (t - t d (x)) — q{x, t)] t>t d {x) 


(25) 


where t d (x) is a variable transport time delay, a ( x, t ) is a wave decay factor, and q (x, t ) is a forcing 
the disturbance w (x, t ) 



dcr 
a {a) 


function due to 


a(x y t) = exp 



u, t - td (x) 4 1 d (a)) da 


q (x, t ) 



iu (o, f - tj, (aQ+ta (o)) ^ 
a (a, t - t d (x) + t d (cr)) 
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The function / (t — t d (x)) represents the wave propagation and must be determined from the boundary condition 
(18) as 

/ (t) = Gu (t) + Ha (L, t) [f {t - t d (L)) - q (L, t)} (26) 


Let T = td (L) be a transport delay period, then Eq. (25) has a series solution 




k~l 


k - 1 


JJ Ha{L,t-mT) 


m— 0 


Gu (t kT) - Y 


Jt Ha (L,t — mT) 


k=0 Lm=0 


g{L,t-kT) (27) 


Similarly, the solution of Eq. (20) is given as 


a i x ) — < a(Ljt) 


0 T < T d (x) 

{t - r d (x)) T > Td ( X ) 


(28) 


where r — tf — t , r d (x) = T — t d (x), and 


g(r) = HPy (Q,t) + HP Y 


k-1 


■ I i-r / 7- . . rr~i\ 

|| Ha i [L, z 4* mi ) 


k—l Lm=Q 


y (0, t + kT) + (F t + HE t ) ll ( t ) 
'k-1 


+ (F t + H E r ) Y 


ft= 1 


|t Ha(L,t + mT) 


,m=0 


+ kT) (29) 


Equations (27) and (29) illustrate the nature of a closed-loop transport system whereby the solutions are expressed 
as time-shifted series of the transport delay period. The resulting DDEs (19) and (22) thus will contain the time-shifted 
series. Therefore, in general a feedback control is difficult to obtain. Nonetheless, there are two simplified solutions 
in the Emit that we need to consider. The first case is associated with a short time horizon when tf < T, for which 
A (0, i) — 0 from Eq, (29) since the system must be causal. Then, we get from Eqs. (27) and (29) 

y (0, t) = Gu (t) - Ha (L> t) q (L t t) (30) 


*{£,*) = -a (Z,t)s(L,t) (31) 

The second case is associated with an infinite time horizon when tf — > oo, for which f (t — T) ~ / (t) and 
g (r - T) ~ g (r). The infinite time horizon case is possible if the time scale of the linear advection equation is much 
greater than the time scale of the ODE. Equivalently, this time scale separation results in jo: (x)| > p (C), where p is 
the spectral radius of the matrix G. Utilizing the series identity 


1 

1 - Ha (L) 


0 


(32) 


we obtain 


y a, t) 


2 / ( 0,0 


Gu (t)-Ha ( L,t) q (L,t ) 
l-Ha(L,t) 


a (L, t) [Gu (t) — Ha ( L , t) q (L, t)] 
1 - Ha (L, t) 


a(L,t)q(L,t) 


(33) 

(34) 


■ (0) A (0, ,) . ' • 1 -p 0 ["» Ml Ml <35, 

1 — Ha \JLtyt) 

In essence, the solution for any t f may be assumed to be bounded between these two cases, we introduce a factor 7 
that represents the effect of the time horizon such that 0 < 7 < 1 , with 7 = 0 corresponding to the first case and 7 — 1 
corresponding to the second case. Then, we axe now able to obtain a feedback control v in terms of the boundary 
control u and the transport states at the boundaries y ( 0 , t) and y (L, t) 


v (t) = — R~ 1 D T Wu — R - 1 D T Va(I- J t) q (L,t) 


( 36 ) 
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where W and V are solutions of the following matrix Riccati equations using a backward sweep method by letting 
jj, = W Au + Va (L, t) q (L, t) 

WC e + CjW-WDR~ I D T W + Q e =0 (37) 

C e T V - WDR- 1 D t V - WS - U = 0 (38) 

with 

C e = C + EG + 7 (F + EH) G a (L, t ) [1 - 7 J?a (£,*)] _1 
Qe = Q + G t PG [1 - 7 Ra (L, t)}~ 2 
S = (F + EH) [1 - -yHa ( L , i?)} -1 
U = G T Pi? [1 - 7 i?a (L, t)] ' ~ 2 

We see that the optimal control solution of a linear advection equation has a form of a Riccati solution. The Riccati 
equation contains modified matrices that incorporate the dynamics of the linear advection equation. The control is a 
state feedback and a disturbance feedforward where the forcing function q ( L , t) is the disturbance that is delayed by 
the time delay T. Thus, the control would not be responsive during this time delay. We note that the control can also 
be written in an output feedback form by noting that 

a (L, t ) q (L, t) = ja (L, t ) y (0, t) - y (L, t) (39) 

so that 

v{t) = — R -1 D t Wu — R~ 1 D T V 7 a (L, t) y (0, t) + R -1 D T Vy (L,t) (40) 

V. Flow Control Application 

We now apply the general theory to a flow control problem to regulate the test section Mach number in a closed- 
circuit wind tunnel. An example is the NASA Ames 1 1 -Foot Transonic Wind Tunnel as shown in Fig. 2. 



Fig. 2 - Closed-Circuit Wind Tunnel 

The fluid transport process in a closed-circuit wind tunnel is a good example of a closed-loop transport process 
whereby the fluid flow is recirculated through a compressor providing the boundary control action. The compressor 
is controlled by two auxiliary dynamical' processes: a drive motor dynamics, and an inlet guide vane (IGV) angle 
dynamics. By controlling the drive motor speed and the IGV angle, the flow in the test section can be set as desired. 

The full nonlinear model of this system is described by Eq. (4) with y (x, t) = \ rh po Tq j and 


A = 


PQC 2 

pA 

(fe-l)T 

pA 


U [l — 


2d 
20 

(fc-l)T 

J 

(fc— 1) 2 Tu 
kp 0 


mu 

2Tq 

P0C 2 U 

To 


U 1 + 


(k~l)T 

To 


B = 


mu$ 

fi _ (fc- 1 !? ] 

2c^ To j 

{fc-l) 2 Tu 3 ^ 

2c 2 


where rh is the mass flow, p is the pressure, T is the temperature , p is the density, u is the flow speed, c is the speed 
of sound, f is a loss factor, A is the cross sectional area, and the subscript 0 denotes the stagnation condition. 
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The compressor provides a boundary control action at x = 0 and x = L that results in a total pressure rise that 
recirculates the flow. Using an empirical model, we obtain the forcing function g in the periodic boundary condition 
(5) for the wind tunnel model as 


m (L, i) 

1+2 YL c ij@ juJ l j ( h ~~ ^2“ — 

where uj is the compressor speed, 0 is the IGV angle, and bi , %•> % are empirical coefficients derived from experi- 
mental compressor performance measurements. 

The drive motors employ a drive motor speed control system that utilizes a variable resistance device known as 
rheostat to control the motor speed. By varying a rotor resistance R r in the rheostat, the motor torque changes. 
The difference between the motor torque and the aerodynamic torque causes the motors to accelerate or decelerate 
according to the following equation 


g (y(L,t); ll) = 


Po (L,t) 



— 


K m RrU; $ (cO s —uj) 


[R s (cu s ~~ lo) + RrUJs] + .Lfwf (cj s - u>y 


- K a [po (o, t) - Po (L, t)} 


(42) 


where J m is the motor inertia, R s is the stator resistance, L s is the stator inductance, w s is the synchronous speed, 
Km is the motor torque constant, and K a is the aerodynamic torque. 

The inlet guide vanes are adjustable with movable trailing edge flaps and are driven by DC field motors that are 
controlled by a field voltage V a according to the following equation 


b + 


K t K e \ x , KiP (A t) « 2 ( L > *) b}c)Cn,e _ K T V a n iP (L, t) u 2 (L, t) b f c}C H 

R a )° ' 2JV 2 "" NR a 2N°- 


(43) 


where b is the viscous friction, K? is the motor torque constant. Kg is the back-emf constant, R a is the shunt 
resistance, N is the gear reduction ratio, ni is the number of inlet guide vanes, and bf, Cf, Cj? t 0, C# are the span, 
chord, hinge moment coefficient derivative, and hinge moment coefficient of the inlet guide vane flaps; respectively. 

Typically, an aircraft model installed inside the test section undergoes a series of angle of attack changes. As the 
aircraft model pitch angle changes, the momentum loss across the test model generates a flow disturbance that travels 
downstream from the test section to the compressor. Because of the fluid transport delay that exists in the system, the 
feedback control at the compressor normally lags the Mach number response in the test section by a time delay since 
the disturbance has not yet reached the compressor before this time. This disturbance causes a pressure perturbation 
that leads to a drop in the test section Mach number. Without any corrective control before the time delay, the Mach 
number will drop below a prescribed tolerance. This motivates us to seek a model predictive control to minimize this 
Mach number deviation by accounting for the time delay. 

The flow perturbation can be modeled by Eq. (17) with y (x, t) = Apo (x, i) as the total pressure perturbation. 
Linearization of Eq. (41) and neglecting the mass flow and total temperature perturbations results in Eq. (17) with the 
following parameters 

a ( x ) — u ( x ) 

^ kM 2 (x)i(x)[l + kM 2 (x)] 

P ^ x,t) ~ 2 [i - m (*)] 


1 + 


k— 1 

1 + A±M 2 ( x) 


> 0 


/ *\ _ fePOiOQ-^oo/ (%) Cd ( [t) ) Am 

{,) 2 A t L m 

where M is the Mach number, Cg is the aircraft model drag coefficient as a function of the pitch angle <p, the overbar 
denotes the nominal condition, the subscript oo denotes test section condition, A m is the model reference wing area, 
L m is the model length. At is the test section area, and / — 1 for x\ < x < X 2 within the test section and / == 0 
otherwise. 
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From the boundary condition (41), we obtain a linearized boundary condition 


V (0, t) = GAu(t) + Hy (L, t) (44) 

where y (0,-f) is the total pressure perturbation at the compressor exit, y (L, t) is the total pressure perturbation at the 
compressor inlet, and Au — Ac o A9 ^ J is the boundary control error vector comprising of the compressor 

speed error Aoj, the IGV flap position error A9, and the compressor speed error integral 0 — Atodr. The error 
integral is designed to ensure a zero steady state error in the compressor speed. In essence, the control scheme is a 
proportional-integral control. The matrix G = | dpo(o,t) q J and If = are the partial derivatives 

evaluated from the boundary condition (4 1), 

The flow disturbance also creates a perturbation in the actuator dynamics of the drive motors and the IGV system, 
thus resulting in the following state equation 

Au - G Au + DAv 4- Ey (0 } t) + F y (L, t) (45) 

where Av = |" AR r AV a j is the augmented control input vector comprising of the augmented drive motor rotor 
resistance ARr, the augmented IGV motor applied held voltage A T / a , and C, D, E, F are partial derivatives evaluated 
from the actuator dynamics, Eqs. (42) and (43), as 



ddj 

Sw 

dib 

m 

0 “ 


dCj 

dR r 

0 


r dc> " 

0po(O,i) 


ddl 

dpo (L>t) 

c = 

M. 

doj 

1 

d& 

dd 

0 

t 

o o 

, D = 

0 

0 

■1 ■ 

. E = 

0 

0 

, E = 

d6 

dvo{L,t) 

0 


Our objective is to design a control law that minimizes the Mach number deviation in the test section. In particular, 
we would like to maintain the Mach number within a required accuracy of ±0.001 at all times. First, we apply the 
following feedback control law to illustrate the problem with time delay 

V (t) = V - R - 1 D T WAu (t) - R.- 1 D T V7a (L, t) y (0, t) + R~ 1 D T Vy (L, t ) (46) 


where v is a nominal control input at the steady state operation of the wind tunnel. 

A control simulation is performed for a test section Mach number Moo ~ 0.6 at a total pressure po,cc = 2116 
psf for a representative test model. The model support on which the test model is mounted has a second-order time 
response as follows 

<t> = <t>d ^1 - e“^coso; m ^ (47) 

where (j>d is the desired pitch angle set point, t m = 5 sec is the time constant, and cu m = 0.1 is the frequency of the 
model support 

To compute the amplification factor a ( x , i) that appears in the control law, we solve the following PDE that is 
completely equivalent to the integral form of Eq. (25) 


1 da 
a(x) dt 


+ +0{x,t )a — 0 


(48) 


subject to an initial condition a (x, 0) — 1 and a boundary condition a (0, t) = 1. 

To solve for Eq. (48), we discretize the wind tunnel into 45 nodes with Ax = 18.074 ft and L — 795 ft We then 
apply the following upwind finite-difference method to compute a (x, t) 


&i } j 4-1 — Q>ij Oii—iAt Ax “b Pi— ljdi—i (49) 

where i = 2, 3, ... , 45 is the index in the x axis, j = 1,2,... , is the time index, and At = 0.01 sec is chosen in order 
to satisfy the Courant-Friedrichs-Lewy (CFL) stability condition 13 


max 


a>iAt 

Ax 


< 1 


The solution surface of a (z, t) is plotted in Fig. 3 
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The dynamic matrices corresponding to the test section Mach number of 0.6 at a compressor speed of 455 rpm are 
computed to be 



~ -0.059709 

1.3223 

0 " 


" -0.32543 

0 


~ -0.0019318 ’ 

c = 

0 

-1.3565 x 10 - 5 

0 

, D= ■ 

0 

7.6358 x 10 “ 5 

, E = 

0 


1 

0 

0 _ 


0 

0 


0 


F .== 


-0.0011615 
3.1673 x 10 “ 9 
0 


G 




-684.50 25.5415 0 


H = 1.6013 


The weighting factors are selected such that P = 0.001, Q — diag (0.01, 0), and R — diag (l, 1 x 10 -7 ). The 
weighting factors Q 22 and R 2 2 for the IGV system are much smaller than the weighting factors Qn and Ru for 
the drive motors since we also want to control the compressor speed as accurately as possible and want the IGV flap 
position to compensate for the total pressure disturbance generated by the test model' 

Fig. 4(a) is a plot of the computed augmented rotor resistance input to the drive motors for 7 •= 0 and 7 = 1 . 
To regulate the test section Mach number, it can be seen that the rotor resistance must decrease as this would result 
in a corresponding increase in the drive motor input torque to compensate for the increase in the total pressure loss 
generated by the test model. Fig. 4(b) is a plot of the computed augmented IGV motor applied field voltage input to 
the IGV system for 7 .== 0 and 7 — L In both Figs. 4(a) and 4(b), the effect of 7 can generally be interpreted as the 
degree of control efforts. Thus, the control effort ranges from a minimum value for 7 — 0 to a maximum value for 
7 — 1. It is also noted that the augmented control inputs exhibit a time delay of about 3 sec, which is the time it takes 
for the disturbance generated by the test model in the test section to propagate downstream to the compressor inlet. 
This time delay is computed to be AT — td ( L ) — id ( ) = 3.56 sec. 

The response of the total pressure perturbation due to the optimal control augmentation is plotted in Fig. 5. As 
can he seen, the total pressure perturbation between the compressor exit at x = 0 and the test section at x — 470 ft 
is effectively controlled to a zero value as desired. A drop in the total pressure occurs immediately right after the test 
model and propagates to the compressor inlet. 
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Fig. 4 - Augmented Control Inputs to Drive Motors (a) and IGV System (b) 



Fig. 5 - Solution Surface of Total Pressure Perturbation 


As noted in Fig. 4, the time delay AT causes the control to be unresponsive during this time delay while the 
aerodynamic flow condition is changing continuously in inherent problem with the Mach 

number feedback control in a wind tunnel To effectively, handle the total pressure disturbance generated by the test 
model, a model predictive control is proposed using the following disturbance feedforward control law 

v (t) = v - R" 1 D t WAu (t + At) - R~ 1 D T Va (£, t + At) q (I, t + At) (51) 

where Au (t + At) and q (L, t -f At) must be computed a priori from a model described by Eq. (17), (44), and (45). 

Thus, the augmented control inputs are evaluated in advance of the error signals and then added to the nominal 
control values. Therefore, the model predictive control is strictly an open-loop control that is based on the computed 
control inputs from the math model using the estimated drag coefficient of the test model and time history of the model 
support. The computed test section Mach number response to the optimal control augmentation using the feedback 
control in Eq. (46) is plotted iri Fig. 6(a) and using the model predictive control in Eq. (51) is plotted in Fig. 6(b). 
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t, sec 


Fig. 6 - Test Section Mach Number Response to Feedback Control (a) and to Model Predictive Control (b) 

It can be seen that the feedback control brings the test section Mach number' closer to its set point, but is unable 
to maintain the test section Mach number to within a required accuracy of ±0.001 at all times. This is due to the 
time delay which causes the aerodynamic flow condition to change without any control input during the first 3.56 sec, 
thereby resulting in the test section Mach number dropping below the tolerance band before the augmented control 
inputs to the drive motors and the IGY system become sufficient to compensate for the total pressure perturbation. In 
contrast, the model predictive control is clearly much more effective in maintaining the test section Mach number well 
within the required accuracy at all times. In all cases, both types of control for 7 == 0 seems to suffer a steady state 
error due to an insufficient control effort. The Mach number response surface is plotted in Fig. 7, showing the Mach 
number distribution throughout the wind tunnel. The test section Mach number response at x == 470 ft is in “trough” 
below the first spike. The two spikes correspond to locations behind the test model and at the compressor inlet. This 
plot illustrates that a boundary control process can usually only control a closed-loop system at either one of the two 
boundaries. In this case, the cost functional is designed to regulate the system response at x ~ 0. 



Fig. 7 - Mach Number Response Distribution 

The compressor speed response £0 the control input augmentation is plotted in Fig. 8 (a) for the feedback control 
and in Fig. 8 (b) for the model predictive control. The feedback control is unable to maintain the compressor speed 
set point to within a required accuracy of ±0.1 rpm. In contrast, the required compressor speed accuracy is easily 
achieved by the model predictive control for 7 = L 
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Fig. 8 - Compressor Speed Response to Feedback Control (a) and to Model Predictive Control To) 

The ICY flap position response to the control input augmentation is plotted in Figs. 9(a) and 9(b) corresponding 
respectively to the feedback control and the model predictive control. The time delay in the IGV flap position response 
is noticeable between the feedback control and the model predictive control. Thus, this suggests that the effectiveness 
of the model predictive control is attributed largely to the control input augmentation to the IGV system. 



Fig. 9 - IGV Flap Position Response to Feedback Control (a) and to Model Predictive Control (b) 

To demonstrate the effectiveness of the model predictive optimal control for 7 — 1 over the enitre subsonic 
operating envelope, we compute the control augmentation for all Mach numbers from 0.4 to 0.9. The results are 
plotted in Fig. 10. As can be seen, the model predictive optimal control is highly effective for all Mach numbers up to 
0.8. The Mach number perturbation is well within the tolerance of ±0.001. However, at a Mach number of 0.9, there is 
a rapid increase in the Mach number perturbation that exceeds the tolerance. The sudden change in the Mach number 
perturbation at a Mach number of 0.9 is directly attributed to the well-known phenomenon of transonic flow where 
many linear perturbation theories break down near a Mach number of unity. The term (3 (x, t) in the linear perturbation 
model reveals that it becomes undefined for M — -1. Thus, the validity of this model excludes the transonic region 
near a Mach number of unity. 
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Fig. 10 - Test Section Subsonic Mach Number Responses 

In general, the accuracy of a model predictive control is predicated upon the goodness of the parameter estimation 
of the disturbance. Any significant deviation from the planned trajectory of the total pressure disturbance will likely 
cause the Mach number to not meet the required accuracy. Thus, to implement a predictive control scheme, it is 
necessary to have a priori knowledge of the time variation of the total pressure disturbance as a function of the drag 
coefficient of the test model. During a pitch polar, the drag coefficient generally varies as a function of the inputs 
parameters that include the pitch angle, the Mach number, and the Reynolds number. A parameter estimation process 
can be established using a recursive least-squares or a neural network algorithm to estimate on-line the drag coefficient 
from the input parameters. In addition, using the knowledge of the time response of the model support system, a 
trajectory of the total pressure disturbance w (t) can be estimated and used to predict the control augmentation from 
the math model. Using the predictive results of the control augmentation, the compressor control would be switched to 
an open-loop mode and begin its actuation simultaneously with the model support system. At the end of the actuation, 
the compressor control would then be switched back to a feedback mode. Because the control is a feedforward scheme, 
stability should not be an issue. The model predictive optimal control thus potentially offers a significant advantage 
over the current feedback control approach which has been demonstrated by simulations and observations to be unable 
to hold the test section Mach number to within a specified accuracy during a continuous pitch polar of the test model. 


VI. Conclusions 

This paper presents some recent results in optimal control of a distributed-parameter system governed by first order, 
quasilinear hyperbolic partial differential equations that is controlled by a boundary control process via a periodic 
boundary condition. The boundary control is further constrained by an ordinary differential equation that models 
an actuator dynamics that exists at the boundary of the system. The resulting coupled hyperbolic partial-ordinary 
differential equation system arises in many physical applications involving transport processes such as traffic or fluid 
flow. Necessary conditions for optimality is derived for this system using the adjoint method which is formulated in 
terms of dual Hamiltonian functions for the partial differential equation and ordinary differential equation systems. 
A linear-quadratic optimal control is developed that results in a two-point boundary value problem involving time- 
shifted solutions due to the periodic boundary condition. By introducing a time horizon parameter, the optimal control 
problem is found to have a Riccati solution. The results are applied to design a model predictive control for the- Mach 
number in a wind tunnel. A feedback control and and a model predictive control are designed to regulate the test 
section Mach number due to a total pressure disturbance generated by a test model pitch motion. The feedback control 
is unable to maintain the test section Mach number to within a prescribed tolerance due to a time delay in the control 
inputs. In contrast, the model predictive control is shown to be highly effective in maintaining the test section Mach 
number to within the required accuracy by relying on a model prediction of the control inputs in advance in order to 
eliminate the time delay effect. 
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